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Summary 

The  first  part  of  this  report  discusses  the  progress  made  in  modeling  the  plume  of  a  high-power 
Hall  thruster  tested  within  the  12V  vacuum  test  facility  located  at  the  Arnold  Engineering  and 
Development  Center  (AEDC),  Tullahoma,  Tennessee.  Over  the  course  of  the  project,  significant 
progress  was  made  in  further  development  of  an  existing  simulation  capability  called  MONACO. 
MONACO  employs  a  combination  of  direct  simulation  Monte  Carlo  and  Particle-ln -Cell 
methods  to  compute  the  plasma  How  in  12V.  New  features  added  to  MONACO  include:  (1) 
analysis  of  solid  baffles  in  12V;  (2)  ability  to  compute  plasma  flow  properties  measured  by  a 
virtual  probe;  and  (3)  implementation  of  a  detailed  electron  fluid  model  to  describe  the  plasma. 
The  new  features  were  assessed  by  direct  comparison  with  data  measured  during  a  Hall  thruster 
test  conducted  in  12V,  and  the  results  were  found  to  be  satisfactory.  The  enhanced  MONACO 
simulation  capability  has  been  transferred  to  AEDC  for  their  future  use. 

The  second  part  of  the  report  describes  progress  made  in  the  simulation  of  a  Hall  thruster 
discharge.  The  need  to  create  an  accurate  model  of  the  Hall  thruster  plasma  generation  was 
identified  as  a  critical  element  for  achieving  accurate  plume  simulations  of  the  thruster  operation 
in  12  V.  The  plasma  discharge  model  simulates  the  unsteady  behavior  of  the  Hall  thruster  using  a 
code  called  H  PH  ALL.  While  H  PH  ALL  has  been  demonstrated  to  produce  physically  accurate 
results  for  time-averaged  properties  such  as  thruster,  there  has  been  little  prior  work  focused  on 
the  unsteady  results  produced  by  the  code.  In  the  work  performed  in  this  project,  detailed 
analysis  was  performed  on  the  unsteady  nature  of  HPHALL  results  with  direct  comparisons  made 
with  unsteady  data  measured  experimentally. 


Part  1--  Modeling  of  Virtual  Diagnostics  for  Hall  Thruster  Plumes  in  a  Vacuum-chamber 


Nomenclature 


e 

electron  charge,  1 .6  x  10"19  C 

k 

= 

Boltzmann  constant  1.38  x  I O’23  J  K'1 

n 

number  density,  m3 

N 

= 

number  of  particles 

m 

— 

mass,  kg 

T 

= 

temperature,  eV 

P 

pressure.  Pa 

specific  impulse,  s 

V 

- 

velocity,  m  s'* 

g 

= 

magnitude  of  relative  velocity,  m  s'1 

&EL 

= 

collision  cross  section,  m2 

A 

“ 

area,  m2 

to 

viscosity  temperature  exponent 

a 

= 

electrical  conductivity,  A  (V  m)'* 

K 

thermal  conductivity,  W  (K  m)"1 

i 

plasma  potential,  V 

2 

= 

electric  field,  V  m'1 

A 

/ 

S= 

current  density,  A  m'~ 

V 

- 

collision  frequency,  s'] 

Q 

= 

ionization  rate  coefficient,  nrf  s'1 

= 

ionization  energy,  eV 

h 

= 

discharge  current,  A 

vd 

discharge  voltage,  V 

V 

electron  velocity  stream  function,  m  2  s 

Subscripts 

a 

= 

anode 

d 

= 

discharge 

e 

electron 

i 

ion 

s 

= 

sampling 

M 

ss 

Maxwellian 

TE 

= 

thruster  exit 

VP 

= 

virtual  probe 

i.i 

Introduction 

Hall  thrusters  are  an  efficient  propulsion  option  for  spacecraft,  with  high  specific  impulses 
making  them  particularly  well  suited  for  low  thrust  missions.  A  primary  concern  regarding  the 
use  of  Hall  thrusters  is  the  effect  of  their  plumes.  Possible  spacecraft  contamination  and 
communications  interference  emphasizes  the  importance  of  plume  analysis  with  regards  to 
spacecraft  integration.  Two  aspects  of  plume  analysis  include  numerical  simulation  and  ground- 
based  experiments  conducted  in  vacuum  chambers.  One  of  these  types  of  facilities  is  the  Arnold 
Engineering  Development  Center  (AEDC)  12V  Vacuum -chamber  (see  Figure  1).  It  is  the  focus 
of  the  present  work  to  develop  simulation  tools  that  can  be  applied  at  AEDC,  specifically  for  the 
analysis  of  Hall  thruster  plumes  in  12V. 


Hall  thruster  plume  modeling  has  been  reviewed  by  Boyd1  where  it  was  determined  that 
hybrid  methods  are  the  most  successful.  In  general,  the  plume  of  a  Hall  thruster  consists  of 
neutrals,  energetic  ions,  and  electrons.  Plume  behavior  is  complicated  by  the  multiple  types  of 
collisions  involved,  such  as  collisions  due  to  thermal  velocity  and  collisions  between  neutrals  and 
ions  with  charge  exchange,  as  well  as  different  physical  phenomena,  such  as  self-consistent 
electric  fields.  Furthermore,  in  ground-based  facilities  the  presence  of  a  background  gas  must  be 
accounted  for.  Computational  methods  are  uniquely  suited  for  plume  analysis  since  different 
physical  models  can  be  interchanged  to  allow  for  varying  degrees  of  fidelity. 

In  this  report,  hybrid  simulation  methods  are  applied  to  an  axisymmetric  domain  in  order  to 
investigate  a  Hall  thruster  plume  in  the  12V  vacuum  chamber,  extending  previous  work  done  on 
the  subject.2,  i  A  direct  simulation  Monte  Carlo  (DSMC)  method4  is  used  to  model  collision 
dynamics,  and  a  Particle-in-Cell  (PIC)  method3  is  used  to  capture  electric  field  effects. 
Additional  capabilities  specific  to  the  needs  of  AEDC  are  also  discussed.  Specifically,  a  virtual 
probe  algorithm  is  implemented  and  compared  against  experimental  data  taken  previously  in  the 
12V  chamber  for  a  high  power  Hall  thruster.  This  algorithm  provides  a  new  capability  in  that 
current  densities  at  discrete  locations  can  be  predicted  more  accurately.  A  detailed  fluid-electron 
model6  is  incorporated  in  place  of  the  standard  Boltzmann  relation  for  modeling  electrons  and  is 
discussed  below.  Comparisons  between  the  fluid-electron  model  and  Boltzmann  relation  model 
are  made  to  determine  effects  on  the  calculation  of  plasma  parameters.  These  comparisons  are 
also  extended  to  determine  the  effect  of  the  electron  model  on  virtual  probe  calculations. 

1.2  Facility  and  Thruster 

The  12V  chamber  is  a  12-ft  diameter  vertically  oriented  vacuum  facility  at  AEDC.  The 
testing  considered  here  was  conducted  in  12V  on  9/13/06  and  9/14/06.  Two  key  features  of  the 
facility  are  the  unique  geometry  and  high  pumping  rates.  The  geometry  is  shown  in  Figure  1 , 
The  floors  and  outer  walls  of  the  chamber  are  liquid  nitrogen  (LN2)  cooled,  held  at  the  cryogenic 
EN2  temperature  of  77K.  T  he  cryo-pumps  located  at  the  top  and  bottom  of  the  chamber  are 
cooled  with  gaseous  helium  (GHe),  held  at  the  cryogenic  GHe  temperature  of  20K.  The  chamber 
has  a  total  pumping  rale  of  about  3-5x1 06  L/s  for  xenon.  This  resulted  in  a  back-pressure  of  1 .4  x 
10"'  torr  for  the  operating  condition  under  consideration  here.  The  current  density  measurements 
utilized  in  this  study  were  obtained  with  Faraday  cup  probes  along  the  chamber  centerline.  For 
further  details,  see  Ref.  7. 

The  thruster  tested  in  this  facility  was  the  Busek  BHT-20K  Hall  thruster.  The  thruster  was 
operated  in  12V  over  a  range  of  conditions  including  a  low  Isp  ,  high  thrust  mode  and  a  high  Isp, 
low  thrust  mode,  the  latter  being  the  focus  of  this  work.  See  Table  I  for  the  specific  operating 
conditions  of  the  high  Isp  mode.  The  BHT-20K  is  a  nominal  20kW  input  power  thruster  designed 
to  produce  1 .0  N  of  thrust  at  2750  s  of  Isp  and  70%  efficiency  at  optimal  operating  conditions,  as 
stated  in  Ref.  8.  For  further  details  of  BHT-20K  operation,  see  the  same  reference. 


Table  I:  BHT-20K  Mall  thruster  operating  conditions 
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Figure  1:  Mesh  for  the  12V  chamber 


1.3  Numerical  Methods 

A.  Overview  and  Collision  Dynamics 

The  numerical  simulations  use  a  hybrid  DSMC-P1C  method  to  model  the  plume.  The  DSMC 
module  handles  collisions.  The  PIC  module  is  used  to  move  the  heavier  particles  that  are 
influenced  by  the  electric  fields  present.  Finally,  the  electrons  are  simulated  using  a  fluid  model, 
incorporating  either  a  simple  Boltzmann  relation  or  a  detailed  model. 

The  DSMC  module  handles  collisions  between  heavier  particles  (Xe,  Xe+,  and  Xe^),  such  as 
neutral-neutral  and  ion-neutral  collisions.  The  DSMC  method  uses  virtual  particles  to  simulate 
collisions  in  rarefied  gas  flows.  The  particles  represent  real  ions  and  neutrals  and  are  grouped  in 
cells  whose  characteristic  lengths  are  shorter  than  a  mean  free  path.  Pairs  of  these  particles  are 
selected  at  random  and  a  collision  probability  is  evaluated  that  is  proportional  to  the  product  of 
the  relative  velocity  and  collision  cross-section.  This  probability  is  compared  to  a  random 
number  to  determine  if  the  collision  occurs.  If  so,  collision  dynamics  are  performed  to  alter  the 
properties  of  the  colliding  particles.  Two  types  of  collision  dynamics  are  relevant  to  Hall  thruster 
plumes:  elastic  (momentum  exchange)  collisions  and  charge  exchange  collisions  (CEX).  Elastic 
collisions  involve  only  exchange  of  momentum  between  participating  particles.  The  l-lall 
thruster  plume  is  confined  to  two  different  types  of  momentum  exchange  collisions,  neutral- 
neutral  collisions  and  neutral-ion  collisions.  For  neutral-neutral  collisions,  the  variable  hard 
sphere  model  is  employed4.  The  cross-section  for  xenon  is: 


aEL(Xe,Xe )  = 


m2  and  =  0.12 


(I) 


2.12X10-1* 

where  g  is  the  relative  velocity  and  co  is  related  to  the  viscosity  temperature  exponent  for  xenon. 
For  neutral-ion  elastic  interactions,  the  cross-sections  measured  by  Miller  et  a I9  are  used: 


oELCXe,Xe*  )  =  (  175.26  —  27.2  log10(p))  x  10  ~20  mz  (2) 

(rEL(Xe,Xe++  )  =  (103.26  -  17.8  log10(tf))  x  1CT20  m2  (3) 

Isotropic  scattering  is  assumed  for  both  types  of  elastic  collisions.  Charge  exchange  collisions 
pertain  to  the  transfer  of  one  or  more  electrons  between  an  atom  and  an  ion.  These  cross  sections 
are  assumed  to  follow  the  same  expressions  for  neutral-ion  elastic  collisions.  However,  it  is  also 
assumed  there  is  no  transfer  of  momentum  accompanying  the  charge  exchange,  since  it  is 
primarily  a  long-range  interaction. 

B.  Plasma  Dynamics 

The  PIC  module  is  used  to  move  the  heavier  ion  particles  that  are  influenced  by  the  electric 

fields,  whereas  the  lighter  electrons  are  modeled  as  a  fluid.  The  PIC  module  determines  the 

charge  density  at  the  nodes  in  the  mesh  based  on  the  proximity  of  each  particle  to  the 
surrounding  nodes.  The  charge  density  is  then  used  to  compute  the  electric  field  at  each  node. 
This  is  accomplished  either  by  incorporating  the  Boltzmann  relation  or  solving  for  the  potential 
directly  using  the  detailed-fluid  model.  The  potential  is  then  differentiated  spatially  to  obtain  the 
electric  fields. 

The  Boltzmann  relation  uses  several  assumptions  applied  to  the  electron  momentum  equation, 
such  as  the  fluid  electron  flow  being  collision  less,  isothermal,  with  no  magnetic  fields  present, 
and  that  the  electron  pressure  obeys  the  ideal  gas  law,  to  arrive  at  the  following10: 

<t>  =  *r.,+  ^T~  ln(^)  (4) 

where  <f>  is  the  plasma  potential,  <j>ref  is  a  reference  potential,  k  is  Boltzmann’s  constant,  e  is  the 
electron  charge,  Tre f  is  the  constant  electron  reference  temperature,  and  ne  and  nref  are  the  local 
electron  number  density  and  reference  electron  number  density,  respectively.  Reference  values 
are  taken  at  the  thruster  exit  plane  and  the  parameters  used  in  this  study  can  be  found  in  Table  1. 
The  assumption  of  quasi-neutrality  is  employed  to  obtain  the  electron  number  density  from  the 
ion  number  densities. 

Most  of  the  strong  assumptions  made  in  deriving  the  Boltzmann  relation  are  questionable 
in  the  plume  under  consideration,  especially  in  the  near-field.  This  is  due  to  strong  gradients  that 
could  make  these  approximations  inaccurate.  Previous  work  on  this  subject  pointed  to  increasing 
the  fidelity  of  physics  models  to  obtain  better  agreement  with  experimental  results'*.  The  detailed 
electron  model  first  proposed  by  Boyd  and  Yinv  increases  the  modeling  fidelity  by  removing 
some  of  the  simplifications  made  with  the  Boltzmann  relation  and  modeling  the  electrons  with 
fluid-type  conservation  equations.  The  relations  are  manipulated  into  useful  forms  for  numerical 


simulation  by  introducing  a  stream -function  (called  the  electron  velocity  potential)  for  the 
electron  continuity  equation  and  assuming  steady  state  for  all  the  equations.  This  results  in  a  set 
of  Laplace-type  equations  with  weak  source  terms,  summarized  as  follows: 


—  n  e  Ci  where  Vip  — 

V(crV0)  =  j(aV:re  +  +  cr^OtO,,))  ■  VT,  +  TeVa  ■  P(te(ne))  +  V<r  ■  )  (6) 

V2T'  =  ■  VT,  +  ^(-j  •  2  +  f  71, GJ  •  rt*r,  +  PtV-^+  3 ^V'intkCT,  -  Tk)  +  (7) 

This  system  is  treated  in  the  current  study  in  the  same  manner  as  in  Refs.  1 1  and  12,  the  results  of 
which  are  summarized  as  follows.  By  treating  the  right-hand  side  as  known  in  Eqs.  (5)-(7),  three 
fundamental  plasma  parameters  are  solved  for,  namely  pe,  and  Te.  This  is  achieved  by 
expressing  the  system  as  a  generalized  Poisson  equation  and  solving  the  system  with  a  finite 
element  solver.  Derivative  calculation  is  handled  by  a  least  squares  method,  also  presented  in 
Ref.  1 1 . 

The  set  of  Eqs.  (5)-(7)  is  used  twice  in  the  current  study.  First,  the  isothermal  detailed  model 
consists  of  solving  Eqs.  (5)  and  (6),  assuming  no  change  in  Te.  Second,  the  conductivity  detailed 
model  consists  in  solving  Eqs.  (5)-(7)  keeping  only  the  first  source  term  in  (7),  in  order  to  ensure 
the  system  of  equations  is  weak.  The  two  detailed  models  increase  the  level  of  accuracy  by 
removing  some  restrictions  that  the  Boltzmann  relation  assumes,  e.g.  currentless  and  isothermal 
electrons. 

C.  Boundary  Conditions 

For  computations  of  a  Hall  thruster  plume  in  12V,  boundary  conditions  must  be  specified  at 
the  thruster  exit  and  along  all  solid  surfaces  in  the  computational  domain,  including  cryo-pumps 
and  baffles,  for  the  DSMC-PIC  method  and  for  the  fluid  electron  model.  First,  the  boundary 
conditions  for  the  DSMC-PIC  method  are  presented,  followed  by  those  for  the  fluid  electron 
model. 

Some  of  the  macroscopic  properties  of  the  plasma  are  required  at  the  thruster  exit,  namely  the 
number  density,  velocity,  and  temperature  of  each  heavy  species  in  the  calculation.  Since  the 
particles  exit  the  thruster  with  an  unknown  radial  velocity  component,  this  results  in  a  velocity 
vector  that  is  not  parallel  to  the  center-line  of  the  chamber.  The  angle  between  the  center-line 
and  velocity  vector  is  referred  to  as  the  divergence  angle,  and  since  it  is  not  known  a  priori,  a 
sensitivity  study  is  performed  that  includes  a  range  of  assumed  values.  Particles  that  exit  the 
thruster  plane  are  assigned  a  radial  velocity  component  that  varies  linearly  with  the  distance  from 
the  center  of  the  thruster  channel.  Therefore,  a  particle  exiting  the  center  of  the  channel  will  have 
no  radial  velocity  component,  while  a  particle  exiting  the  channel  near  its  inner  or  outer  wall  will 
have  the  highest  radial  velocity  component  and  the  angle  of  the  velocity  vector  with  respect  to 
the  center-line  for  that  particle  will  be  the  divergence  angle  specified. 

The  other  thruster  exit  macroscopic  properties  are  determined  in  general  via  a  combination  of 
analysis  and  estimation  in  order  to  match  experimental  operating  conditions  in  the  same  manner 
as  Ref.  3.  The  resulting  DSMC-PIC  input  parameters  produce  thruster  exit  conditions  matching 
the  experimental  conditions  to  within  5%.  This  has  been  noted  in  the  same  work  as  a  source  of 
uncertainty  in  comparison  with  experimental  data.  The  experimental  operating  condition  under 
consideration  here  is  a  high  Isp,  high  voltage  configuration,  as  shown  in  Table  1.  See  Table  2  for 


the  listing  of  thruster  exit  simulation  parameters  used  to  replicate  those  operating  conditions.  It 
should  be  noted  that  the  BHT-20K.  has  undergone  almost  no  detailed  characterization  at  the 
thruster  exit  and  so  the  determination  of  appropriate  thruster  exit  conditions  is  much  more 
difficult  than  for  the  Aerojet  BPT-4000  thruster  studied  in  12V  previously2. 

The  DSMC-PIC  boundary  conditions  pertaining  to  the  Boltzmann  model  and  the  detailed- 
fluid  model  differ  with  regard  to  the  ion  number  density  and  velocity  in  order  to  obtain  far-field 
agreement  between  the  models,  as  characterized  by  the  axial  velocity  along  the  center-line.  The 
flow  velocities  at  the  exit  are  different  in  order  to  obtain  the  far-field  agreement  and  the  number 
densities  are  adjusted  to  maintain  the  desired  mass  flow  rate.  This  agreement  is  obtained  as 
shown  in  Figure  2,  in  which  the  thruster  exit  is  located  at  z  =  3  m. 

Concerning  the  different  types  of  solid  surfaces  in  12  V,  each  one  is  assumed  to  have  a  plasma 
potential  of  zero,  including  the  baffles,  which  are  electrically  grounded.  All  ions  that  collide 
with  any  wall  are  neutralized.  When  particles  strike  the  cryo-pump  surfaces,  a  fraction  of  them 
are  pumped  away.  This  process  is  characterized  by  a  sticking  coefficient  (a  value  of  0.8  is  used 
for  the  present  study).  For  the  particles  scattered  back  into  the  flow  field  from  all  surfaces, 
diffuse  reflection  is  assumed  which  is  characterized  by  the  surface  temperature. 

The  facility  back-pressure  is  modeled  through  static  background  particles.  Each  cell  contains  a 
few  particles  with  velocities  sampled  from  a  zero-centered  Maxwellian  velocity  distribution 
function.  These  particles  participate  in  collisions  with  plume  particles  and  change  the  velocities 
of  other  particles,  but  their  positions  and  velocities  do  not  change.  The  back-pressure  value  for 
these  simulations  is  set  to  the  value  recorded  during  the  experiment,  1 .4  x  1G'7  torr. 


Species 

Number  Density  (  1  / 

■ ii 

Axial  Velocity  (  m  /  s  ) 

Temperature  (  K ) 

Boltzmann 

Model 

Detailed 

Model 

Boltzmann 

Model 

Detailed 

Model 

Boltzmann 

Model 

Detailed 

Model 

Xe 

1 .28  x  10IK 

1.28  x 
10,s 

281.3 

281.3 

750 

750 

Xe+ 

1.18  x  1017 

1.39  x 
1017 

27792.0 

25792.0 

23188 

23188 

Xe~ 

2.96  x  1016 

3.42  x 

1016 

39304.0 

36475.0 

23188 

23188 

Table  2:  DSMC-PIC  input  parameters  at  the  thruster  exit 


Figure  2:  Profiles  of  Xe+  axial  velocity  along  the  plume  centerline 

The  boundary  conditions  for  the  detailed-fluid  model  must  be  determined  in  order  to  solve 
F,qs.  (5)  -  (7)  above.  Since  each  equation  is  Laplace-like,  each  one  requires  either  a  Dirichlet 
(direct)  value  specified,  or  a  von  Neumann  (gradient)  value  specified.  The  potential,  $  the 
electron  temperature,  Te,  and  the  quantity  neve  are  specified  as  either  direct  or  gradient  at  each 
boundary  in  the  simulation.  The  value  to  which  each  is  set  is  specified  in  Table  3.  Each 
boundary  condition  is  direct  except  for  the  following:  1)  each  boundary  condition  for  the  axis  of 
symmetry  is  a  gradient-type  condition,  and  2)  the  wall  of  the  thruster  has  a  gradient-type 
condition  for  the  electron  temperature.  Due  to  the  grounding  of  all  the  surfaces  in  the  chamber, 
the  plasma  potential  condition  is  direct  and  set  to  zero  for  the  chamber  walls  and  set  to  reference 
values  for  the  thruster  exit  plane  and  cathode.  The  electron  temperature  is  set  to  nominal  values 
based  on  previous  work  and  numerical  stability.  The  boundary  conditions  for  neve are  determined 
using  conservation  of  current  as  follows: 


Jd  -  h  +  h  (8) 

T  —  h  *  (9) 

u 

JiLjtatkQ&M  —  .  (10) 

where  (8)  and  (9)  are  used  to  compute  current  density  at  the  thruster  exit  and  (  !0)  computes 
current  density  at  the  cathode,  from  which  neve  is  found  by  dividing  by  the  elementary  charge  e. 
The  sign  convention  results  from  electrons  flowing  out  of  the  cathode,  towards  the  anode. 
Therefore,  the  electron  velocity  stream  function  must  be  positive  at  the  cathode  and  negative  at 
the  thruster  exit. 


Boundary 

Condition 

Variable 

Thruster 

Cathode 

Symmetry- 

line 

(all  gradient- 

type) 

Chamber 

Walls 

Exit  Plane 

Body 

dW 

-1.79  x  !021 

m'V 

0 

7.32  x  1022 
m’V 

0 

0 

4 

50  V 

0  V 

5  V 

0  V  nf1 

0  V 

Te 

10  eV 

0  eV  m'1 
(gradient) 

2  eV 

OeV  m'1 

1  eV 

Table  3:  Detailed  fluid-electron  boundary  conditions 
D.  Virtual  Probe  Modeling 

A  virtual  diagnostics  algorithm  is  implemented  in  order  to  compare  current  density  predictions 
to  measured  data.  The  virtual  probe  is  analogous  to  a  Faraday  cup  probe,  whereby  ions  crossing 
the  virtual  probe  surface  are  counted.  This  flux  is  then  converted  to  current  density. 

Once  the  number  of  ions  that  pass  through  the  virtual  probe  face  over  the  course  of  the 
simulation  is  ascertained,  this  can  be  converted  to  a  current  through  the  following  formula: 

/  =  (11) 

where  I  is  the  calculated  ion  current,  qu  is  the  charge  of  particle  k,  and  ts  is  the  total  sampling 
time.  Once  the  current  is  computed,  it  is  divided  by  the  area  of  the  virtual  probe  collecting 
surface  to  calculate  a  current  density.  This  allows  for  a  direct  comparison  to  experimental  results 
while  still  employing  an  axisymmetric  domain.  'This  method  is  versatile  in  that  it  can  account  for 
a  variety  of  probe  locations  and  orientations,  as  well  as  probe  operating  parameters  such  as 
potential  biasing,  without  affecting  the  flow-field  solution. 

Preliminary  assessment  of  the  algorithm  is  performed  for  a  virtual  probe  with  no  potential 
bias,  aligned  perpendicular  to  the  thruster’s  acceleration  channel  at  a  location  10  cm  downstream 
of  the  thruster  exit  plane.  Here,  the  current  density  can  be  effectively  approximated  with  the 
following  macroscopic  formula: 


;=£i<nu<r)i  (i2) 

where  u  is  the  average  axial  velocity  component  of  the  ions,  q\  is  the  elementary  charge  of  the 
ion,  and  n  is  the  local  number  density  of  ions.  The  radial  component  of  the  ion  velocity  can  be 
neglected  since  it  is  much  smaller  than  the  axial  component.  The  approximated  current  density 
using  Eq.  (12)  agrees  with  the  prediction  using  the  virtual  probe  to  within  1 .2%. 


1.4  Plume  Simulation  Results 


The  simulations  presented  use  a  total  of  250.000  particles  at  steady  state  over  a  domain  of 
1 853  triangular  cells.  The  simulation  runs  for  800,000  time-steps  to  reach  a  steady  state  and  then 
for  another  500,000  time-steps  to  sample  macroscopic  data.  The  time-step  size  is  2  x  10"(> 
seconds,  resulting  in  a  total  sampling  time  of  1  second.  Cases  are  run  at  divergence  angles  of 
10°,  20°,  and  30°, 

A.  Flowfield  Comparisons 

Shown  below  in  Figures  3-8  are  field  data  resulting  from  the  simulations  described  above. 
Qualitative  effects  of  divergence  angle  are  observable  on  each  figure,  with  10°  and  30°  cases 
shown  side  by  side.  The  difference  in  current  density  predictions  can  be  seen  in  Figures  3  and  5. 
Mote  that  the  current  density  computed  with  the  conductivity  model  is  not  included  here  due  to 
its  similarity  with  the  current  density  computed  by  the  isothermal  model.  The  effect  of  using  the 
detailed  model  is  shown  to  distribute  the  current  density  throughout  a  greater  part  of  the  domain 
when  compared  to  the  Boltzmann  model.  This  effect  is  due  to  the  mathematical  nature  of  the 
Boltzmann  model  in  that  it  restricts  the  dynamic  range  of  the  potential.  This  can  be  seen  directly 
through  the  plasma  potential  (Figures  4,  6,  8),  where  the  Boltzmann  model  restricts  the  dynamic 
range  in  potential  to  around  25  V,  whereas  the  detailed  models  allow  for  variation  of  around  50 
V.  Finally,  Figure  7  shows  the  qualitative  difference  between  the  Boltzmann/isothermal  models 
and  the  conductivity  model:  the  electron  temperature  is  held  at  a  constant  2  eV  in  the  former 
models  and  solved  for  directly  in  the  latter  model.  This  effect  can  be  seen  in  comparing  Figures 
6  and  8:  the  conductivity  model  predicts  even  sharper  gradients  in  the  plasma  potential  than  the 
isothermal  model. 

The  effect  of  divergence  angle  choice  can  be  observed  most  explicitly  in  the  Boltzmann  model 
results,  and,  to  a  lesser  degree,  in  the  results  of  the  detailed  models.  A  smaller  value  of 
divergence  angle  results  in  a  more  beam-like  plume.  This  behavior,  seen  in  Figures  3  and  5, 
indicates  the  predicted  current  densities  are  highly  dependent  upon  the  choice  of  divergence 
angle.  Since  divergence  angle  is  strictly  a  numerical  feature  not  known  for  the  BHT-2QK, 
analyzing  a  range  of  angles  allows  us  to  fine-tune  the  results. 
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Figure  3:  Current  density  using  the 
Boltzmann  fluid-electron  model 


Figure  4:  Plasma  potential  using  the 
Boltzmann  fluid-electron  model 
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Figure  5 :  Current  density  using  the 
isothermal  detailed  fluid-electron  model 


Figure  6:  Plasma  potential  using  the 
isothermal  detailed  fluid-electron  mode! 
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Figure  7;  Electron  temperature  comparisons  -  Figure  8;  Plasma  potential  using  the 
conductivity  and  Boltzmann/isothermal  Models  conductivity  detailed  fluid -electron  model 


B.  Experimental  Comparisons 

Shown  below  in  Figures  9-11  are  comparisons  between  experimental  measurements  and 
various  simulations.  Measured  data  are  reported  in  Ref.  13  and  data  acquisition  methods  are 
described  in  Ref.  8  corresponding  to  the  experiments  performed  in  12V  for  the  Busek  BHT-20K 
Hall  thruster  at  the  operating  condition  under  consideration  here.  Model  assessment  is  performed 
through  comparison  with  data  measured  by  a  set  of  Faraday  cup  probes  whose  locations  are 
shown  in  Figure  1.  Note  that  two  probes  are  used  with  two  different  thruster  locations,  resulting 
in  four  points  of  experimental  measurements,  i.e.  probes  1  &  2  and  1’  &  2’.  A  third  probe  was 
present  upstream  of  probes  1  and  2  but  not  shown  here  due  to  its  shorting  out  at  the  beginning  of 
the  test.  General  trends  seen  in  Figures  9-11  capture  the  tendencies  observed  in  experimental 
measurements;  note  that  discrepancies  are  tabulated  in  Table  4  for  the  30°  divergence  angle  case 
using  the  detailed  models  only:  these  models  provide  the  closest  comparison  when  considering 
the  discrepancies  as  a  percent  difference  between  experimental  measurements  and  computational 
predictions. 

The  error  bars  shown  are  based  on  probe  shielding  interference.  Probe  1  was  shielded  by  the 
shorted  probe  not  shown,  whereas  probe  2  was  shielded  by  both  the  shorted  probe  and  probe  I . 
Assuming  a  maximum  shielding  based  on  geometry,  the  shorted  probe  shields  a  maximum  of 
14%  of  the  ion  current  for  probe  1,  while  the  shorted  probe  and  probe  I  shield  a  maximum  of 
37%  of  the  ion  current  for  probe  2. 


The  two  prediction  methods  utilized  are  the  virtual  probe  algorithm  described  above,  and 
evaluation  based  on  the  macroscopic  properties  of  the  plasma.  The  macroscopic  properties  are 
determined  by  taking  averages  of  particles  within  a  cell,  whereas  the  virtual  probes  are 
implemented  only  at  the  experimental  measurement  locations;  as  such,  the  virtual  probe 
predictions  are  presented  below  as  discrete  points.  The  differences  between  the  macroscopic 
model  and  the  virtual  probe  will  be  discussed  below.  The  virtual  probe  data  are  color  coded  to 
the  fluid-electron  model  utilized  for  that  case. 
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Figure  9:  Current  density  comparisons,  10°  divergence  angle 
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Figure  10:  Current  density  comparisons,  20°  divergence  angle 


Figure  II;  Current  density  comparisons,  30°  divergence  angle 
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Table  4:  Tabulated  current  densities— 30°  divergence  cases 


When  comparing  the  predictions  to  the  experimental  measurements,  agreement  depends 
highly  upon  mode!  choice  and  divergence  angle.  The  Boltzmann  model  virtual  probe — 10° 
divergence  angle  case  results  in  a  prediction  which  only  captures  qualitative  trends  in  the  data, 
over-predicting  measured  current  density  by  almost  a  factor  of  JO,  whereas  using  either  of  the 
detailed  models  for  a  30°  divergence  angle  results  in  discrepancies  of  0,48-40,1%,  Most  of  these 
discrepancies  fall  within  the  uncertainty  shown  for  probe  interference. 

C.  Comparing  the  Virtual  Probe  and  Maxwellian  Models 

The  significant  difference  between  the  predicted  current  densities  using  the  macroscopic 
assumption  or  the  virtual  probe  model,  shown  in  Figs.  9-11,  can  be  explained  by  analyzing  the 
velocity  distributions  that  each  model  predicts.  The  Maxwellian  velocity  distribution  function 
(VDF)  in  the  axial  direction  is14: 

/*  < Cl > dC,  =  «p (“ci)  dCt  (13) 

The  VDF  predicted  by  the  virtual  probe  model  is  determined  by  computing  a  histogram  based 
on  the  velocities  of  the  recorded  particles  passing  through  the  probe  face.  A  current  distribution 
function  (CDF,  G)  is  also  computed  for  a  closer  representation  of  the  difference  in  models.  The 
CDF  assuming  a  Maxwellian  VDF  is  computed  by  Fq.  (14).  The  CDF  is  computed  for  the 
virtual  probe  using  Fq.  (15).  Shown  below  in  Figures  12  and  13  are  the  VDF  and  CDF  predicted 
by  each  model.  By  integrating  the  CDF,  the  predicted  current  density  is  obtained,  therefore  the 
total  area  under  each  CDF  curve  represents  the  current  density  predicted  for  each  method.  A 
visual  comparison  of  each  VDF  shows  that  the  velocity  distribution  corresponding  to  the  ions 
crossing  the  virtual  probe  plane  is  non-Maxwell  tan:  the  beam  ions  are  tightly  centered  around  29 
km/s,  with  a  small  group  of  CEX  ions  centered  at  4  km/s.  The  corresponding  CDF  for  each 
method  further  exhibits  the  inaccuracy  of  the  Maxwellian  assumption  at  the  probe  locations:  the 


integrated  current  density  based  on  a  Maxwellian  VDF  is  much  greater  than  the  integrated 
current  density  based  on  the  virtual  probe.  Therefore,  a  VDF  based  on  the  equilibrium 
assumption  is  inappropriate.  Since  the  virtual  probe  is  based  on  discrete  particle  detection,  it 
captures  non-equi librium  effects  directly.  The  result  is  a  current  density  prediction  that  falls 
within  the  uncertainty  associated  with  the  experimental  set-up. 


Gx  (Ct)—  Jiiuni  CiA  f  Ct  )dCt 
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Figure  t2:  Maxweliian  and  virtual  probe  velocity 
distribution  functions  at  probe  divergence  angle 
of  10°,  conductivity  model  case 


Figure  13:  Maxwellian  and  virtual  probe  current 
distribution  functions  at  probe  I,  divergence  angle 
of  10°,  conductivity  model  case 


1.5  Conclusions 

A  genera]  purpose,  hybrid  DSMC-P1C  code  has  been  enhanced  with  a  new  simulation 
tool  for  predicting  current  densities  and  a  new  fluid-electron  model  which  enhances  physical 
modeling  fidelity.  The  new  virtual  probe  model  provides  enhanced  accuracy  in  comparison  to 
experimental  measurements  and  is  flexible  enough  to  provide  predictions  at  locations  that  were 
previously  lacking  instrumentation.  This  new  level  of  modeling  has  increased  the  understanding 
of  the  test  case  under  consideration.  The  inclusion  of  more  robust  fluid-electron  models 
alleviated  restrictions  that  were  noted  in  previous  work3  and  provided  higher  correlation  with 
experiment.  The  result  of  these  new  modeling  applications  showed  excellent  agreement  to 
experimental  measurements,  improving  prediction  capabilities  from  merely  qualitative  and 
deviations  from  experimental  data  as  large  as  a  factor  of  1 0,  to  predictions  which  fall  within 
uncertainties  of  the  measurements. 

Further  refinement  of  thruster  exit  plane  boundary  conditions  and  augmentation  of  the  virtual 
probe  model  to  account  for  probe  operating  parameters,  such  as  bias  voltage,  are  areas  of  future 
study  to  better  understand  the  physical  situation  in  the  plume.  Other  areas  of  study  include 
refinement  of  the  detailed  fluid-electron  model  in  order  to  incorporate  all  the  source  terms 
associated  with  Eq.  (7)  while  maintaining  numerical  stability.  Specifically,  the  Joule  heating 
term  could  play  an  important  role  in  capturing  higher  order  effects. 
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Part  2—  Spectral  analysis  of  simulated  Hall  thruster  discharge  current  oscillations 

Nomenclature 

I  d  =  discharge  current 

F  th  =  force  due  to  thrust 

f  =  frequency 

At  =  base  (ion)  time  step 

At  e  =  electron  subcycle  time  step 

m  a  =  mass  of  species  a 

n  a  =  number  density  of  species  a 

q  =  elementary  charge 

v  a  =  collision  frequency  for  phenomenon  a 

p  ez  =  axial,  cross-field  electron  mobility 

a  -  location-based  scaling  coefficient  for  electron  mobility  modeling 
to  a  =  plasma  frequency  of  species  a 
Z  =  charge  number 
B  i  =  magnetic  field,  i  component 

2.1  Introduction 

Hall  thrusters  are  known  to  generate  oscillation  modes  as  a  result  of  physical  processes, 
some  of  which  are  known  and  others  which  are  still  being  observed  and  studied  {Choueiri, 
2001 },  making  time-varying  modeling  an  important  research  topic.  For  the  past  decade,  the  Hall 
thruster  model  H  PH  ALL  has  been  the  focus  of  much  research  and  development.  H  PH  ALL 
performs  an  axisymmetric  simulation  of  the  plasma  within  the  thruster  discharge  chamber  and 
near-field  plume  which  employs  both  fluid  and  particle-in-cell  (PIC)  numerical  methods  {Fife, 
1998}.  The  code  has  been  found  to  be  potentially  effective  in  creating  either  time-averaged 
outputs  of  performance  data  which  could  be  used,  for  example,  as  source  terms  into  an  erosion 
mode  {Hofer  et  al.,  2007},  or  to  observe  and  analyze  the  time-varying  nature  of  the  detailed 
evolution  of  the  internal  plasma  of  a  Hall  thruster  at  small  timescales  {Parra  et  al.,  2006}. 
However,  it  seems  that  the  time-varying  results  have  undergone  less  scrutiny  due  to  the  fact  that 
while  time-averaged  results  are  tunable  for  accurate  steady-state  comparisons  to  experiment, 
time- varying  observations  can  be  plagued  with  inconsistencies.  It  is  this  time-varying  nature  and 
subsequent  phenomena  of  time-varying  operation  which  is  the  focus  of  this  report. 

This  report  provides  a  brief  introduction  to  the  Hall  thruster  under  consideration  and  the 
simulation  code.  The  methods  of  comparison  for  the  present  work  are  outlined  whereby 
simulated  results  of  both  performance  and  spectral  analyses  are  compared  to  experimental  results 
for  assessment  purposes.  Assessment  of  the  code's  ability  to  replicate  accurate  performance 
output  is  shown  to  be  successful  for  time-averaged  data  for  the  Hall  thruster's  operating 
conditions  of  interest.  However,  some  of  the  numerical  parameters  used  in  the  simulation  of 
these  operating  conditions  are  shown  to  have  adverse  effects  on  the  response  of  physical  models. 
Varying  responses  due  to  physical  models  are  then  shown  to  have  an  effect  on  time- varying 
characteristics  of  the  performance  output.  The  results  and  disparities  of  this  assessment 
demonstrate  the  consequent  need  for  a  more  detailed  look  at  numerical  schemes  of  the  physical 


models  as  well  as  the  numerical  parameters  employed.  Steps  are  taken  to  vary  simulation  particle 
population,  time  step,  and  the  length  scales  of  certain  physical  models  in  order  to  study  the  way 
in  which  time-varying  data  is  affected. 

2.2  Technical  Approach 

A.  Hall  thruster  and  Experiments 

The  experiments  used  for  the  comparison  purposes  of  this  study  were  performed  by  Reid 
using  a  6  kW  laboratory  model  Hall  thruster  {Reid,  2009}.  Experiments  were  performed  in  the 
Large  Vacuum  Test  Facility  (LVTF),  a  6  m  diameter  by  9  m  long  cylindrical,  stainless  steel 
chamber  with  seven  cryopumps,  at  the  University  of  Michigan’s  Plasm adynamics  and  Electric 
Propulsion  Laboratory  (PEPL).  The  thruster  was  operated  on  a  100  kW  power  supply  using  a 
separate  1  kHz  RC  discharge  filler  for  protection.  Commercially  available  power  supplies  were 
used  to  power  the  cathode  heater,  cathode  keeper,  and  magnet  circuitry.  Research  grade  xenon 
propellant  (99.999%  pure)  was  used  for  cathode  and  anode  supply. 

In  order  to  gather  time-averaged  discharge  current  data,  a  calibrated  current  shunt  and 
multimeter  were  used  in  correlation  with  the  main  discharge  power  supply  read-out.  Time- 
varying  discharge  current  data  was  gathered  using  a  commercially  available,  high-speed  current 
shunt  rated  for  100  kHz  and  placed  on  the  cathode  return  line  in  between  the  thruster  and  RC 
filter. 

B,  Numerical  Model 

The  computer  code  HPHALL  performs  an  axisymmetric  simulation  which  is  commoniy 
referred  to  as  "hybrid- PIC",  utilizing  fluid  approximation  equations  in  the  treatment  of  electrons 
and  a  particle-in-cell  (PIC)  method  in  the  treatment  of  heavy  species,  namely  singly-  and  doubly- 
charged  xenon  ions  and  atoms.  The  electron  equations  are  solved  at  a  smaller  time  step,  called 
the  electron  subcycle,  in  order  to  simulate  how  the  speed  of  electrons  is  much  faster  than  that  of  a 
heavy  species.  This  allows  for  electron  fluid  equations  to  be  fully  solved  and  settled  in  between 
xenon  particle  updates.  Hybrid  fluid-particle  methods  have  been  shown  to  be  successful  in  Hall 
thruster  plume  studies  {Boyd,  2001}  and  are  computationally  cheaper  than  fully  kinetic  methods 
{Szabo,  2001 }.  The  code  was  originally  created  by  Fife  and  Martinez-Sanchez  {Fife,  1998}  and 
then  later  advanced  by  Gamero-Castano  and  Katz  to  include  such  upgrades  as  a  more  detailed 
sheath  model  and  a  sputtering  yield  algorithm  for  the  use  as  an  erosion  mode!  {Gamero  et  al., 
2005}.  Parra  et  al.  also  made  certain  advances  by  way  of  improving  such  algorithms  as  heavy 
particle  modeling,  electron  mobility  modeling,  plasma  weighting,  and  ionization  models,  among 
others  {Parra  et  al.,  2006a,  2006b}.  Further  development  and  corrections  were  performed  by 
Hofer  et  al.  on  the  heavy  particle  modeling,  erosion  sub-model,  and  electron  mobility  physics, 
continuing  the  development  of  the  code  to  the  present  version.  {Hofer  et  al.,  2006,  Hofer  et  al., 
2007a,  2007b,  Hofer  et  al.,  2008}  Through  this  past  work,  the  code  has  a  favorable  history  of 
presenting  good  agreement  with  macroscopic  properties,  such  as  discharge  current  and  thrust,  as 
well  as  local  properties,  such  as  plasma  density,  plasma  potential,  and  electron  temperature.  An 
example  of  this  time-averaged  output  of  the  code  can  be  seen  in  Figure  I  which  is  a 
representation  of  electron  number  density  at  optimal  operating  conditions,  detailed  in  the  next 


section.  For  this  type  of  two-dimensional  contour  plot,  the  anode  is  located  on  the  left  side  and 
the  near-field  plume  is  located  on  the  right  side  of  the  domain. 
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Figure  1 :  Axisymmetric,  time-averaged  electron  number  density  contours 


3.3  Results 

A.  Code  Assessment 

In  order  to  assess  the  code  using  experimental  results,  simulations  are  performed  at 
several  different,  physical  operating  conditions  which  are  nominal  to  the  thruster's  operation  and 
explained  in  detail  later.  Simulations  are  performed  by  injecting  neutral  xenon  into  the  domain, 
with  plasma  physics  turned  off,  for  approximately  1  ms  of  simulated  time  followed  by 
approximately  4  ms  of  simulated  time  where  the  plasma  physics  are  turned  on.  For  example, 
when  operating  at  the  time  step  At  =  5e-8  s,  20,000  iterations  of  neutral  injection  are  followed  by 
80,000  iterations  of  plasma  simulation.  At  these  operating  conditions,  and  using  approximately 
300,000  particles,  this  simulation  requires  about  12  hours  of  computer  run  time.  All  comparisons 
are  performed  with  electron  subcycles  of  either  0.01  At  or  0.00 1  At  which  have  little  effect  on 
either  the  time-averaged  or  time-varying  solutions  of  the  simulation.  The  electron  subcycle  is 
discussed  further  in  the  time  step  study  section. 


Flowrate 

MA) 

Experiment 

U  (A) 

HPHALL 

%  Diff. 

Fth  (N) 

Experiment 

Fth  (N) 
HPHALL 

%  Diff. 

10  '"Vsec 

9.06 

8.82 

2.6 

0.185 

0.175 

5.4 

20  mVs*c 

19.94 

19.80 

0.7 

0.398 

0.394 

1.0 

20 

19.94 

19.88* 

0.3 

0.398 

0.386* 

3.0 

30  "‘Vsec 

33.80 

32.14 

4.9 

0.620 

0.599 

3.4 

Table  I:  Hall  thruster  operating  conditions  for  both  experiment  and  simulation  at 
a  discharge  voltage  of  300  V,  *  indicates  facility  background  pressure  included. 


Nominal  experimental  operating  conditions  include  discharge  voltages  of  150  V,  300  V, 
and  600  V  at  mass  flow  rates  of  10,  20,  and  30  mg/s.  Simulations  are  run  at  all  operating 
conditions  but  are  focused  on  the  optimal  operating  conditions  of  300V,  20  mg/s.  Simulated, 
time-averaged  performance  data  are  compared  to  experimental  performance  data.  The  two  major 
comparisons  being  discharge  current,  ld,  and  thrust,  Ft)„  and  show  that  simulated  results  are  very 
accurate.  This  comparison  of  varying  flowrate  for  the  optimal  discharge  voltage  of  300  V  can  be 
seen  in  Table  1.  Included  in  Table  1  is  the  condition  of  added  facility  background  pressure,  an 
algorithm  which  injects  neutral  particles  at  the  simulation  boundary.  This  method  brings 
performance  output  closer  to  experimental  results  and  ultimately  simulates  a  more  accurate 
portrayal  of  the  real  system.  The  reason  for  the  best  agreement  at  optimal  conditions  is  because  of 
the  internal  tuning  of  electron  mobility  physics  within  the  simulation  which  has  been  pre-set  to 
show  the  best  agreement  with  these  optimal  conditions,  the  conditions  that  the  Hall  thruster  is 
designed  for  and  at  which  it  would  most  likely  be  operated.  This  time-averaged  comparison 
shows  very  good  agreement  between  simulation  results  and  experimental  measurement  at  the 
operating  conditions  of  interest. 

Discharge  oscillations  from  both  simulation  and  experiment  are  also  compared  via  power 
spectra.  Simulated  trends  can  be  seen  in  the  Hall  thruster  discharge  current  spectral  power 
analysis  found  in  Figure  2  and  show  relatively  good  agreement  with  experimental  data  {Reid, 
2009},  also  found  in  Figure  2  (using  the  same  color  scheme  for  the  purpose  of  comparison), 
though  the  agreement  is  far  from  exact.  The  noisy  simulated  data  has  been  processed  through  a 
basic,  one-dimensional  smoothing  algorithm  for  the  purely  observational  purposes  of  comparing 
relative  magnitudes  and  trends.  In  terms  of  absolute  values,  it  is  found  that  the  model  simulates  a 
much  larger  magnitude  discharge  current  oscillation  than  the  experimental  data.  At  optimal 
operating  conditions,  data  from  Reid  {Reid,  2009}  report  a  standard  deviation  of  the  discharge 
oscillations  as  approximately  8%  of  the  mean  discharge  current  whereas  the  simulations  show  an 
approximately  8-25%  oscillation  (depending  on  time  step  size  and  mobility  model,  as  will  be 
shown  later),  which  translates  to  approximately  25-50%  in  peak-lo-peak  terms.  The  familiar 
breathing  mode  oscillation  appears  in  the  range  of  10-30  kHz  in  both  spectra.  A  higher 
frequency  mode  is  also  clearly  present  at  around  70-100  kHz  although  its  source  is  not  clear. 
While  the  fact  that  this  mode  is  seen  in  both  data  sets  seems  like  a  potentially  important  finding, 
it  will  be  shown  that  the  computational  result  is  sensitive  to  physical  modeling  and  numerical 
parameters.  These  results  give  a  glimpse  into  the  details  of  such  time-varying  phenomena  which 
will  affect  the  validity  of  time-varying  comparisons. 
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Figure  2:  Power  spectral  density,  simulation  (top)  and  experiment  (bottom). 


B.  Effects  of  Electron  Mobility  Modeling 


The  first  phenomena  which  must  be  discussed  are  the  electron  mobility  models  currently 
available  within  H  PI  I  ALL  and  their  effects  on  time-varying  solutions.  Cross-field  electron 
mobility  is  modeled  in  the  code  as  {Hofer  et  ah,  2008}: 


i/eme 


(i) 


Vez  - 


a  measure  of  the  propensity  for  electrons  to  cross  magnetic  field  lines  axially  rather  than  stay  in 
an  azimuthal  drift  within  the  discharge  channel.  The  effective  electron  collision  frequency,  14,  is 
given  by 
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with  ven  as  the  electron-neutral  collision  frequency,  vei  as  the  electron -ion  collision  frequency, 
vy-aii  as  the  collision  frequency  of  the  electrons  with  the  walls,  and  Hiohm  as  a  collision  frequency 
due  to  anomalous  Bohm  diffusion  for  electrons.  The  anomalous  Bohm  diffusion  parameter  is  a 
user-augmented  parameter  which  holds  the  purpose  of  matching  simulated,  effective  collision 
frequencies  with  those  of  experimental  measurements.  Anomalous  Bohm  diffusion  is  modeled  in 
the  code  as 
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where  a  is  an  arbitrary  parameter  adjusted  to  match  experimental  results  by  providing  for  the 
necessary  location-based  cross-field  diffusion  values.  For  example,  the  case  for  classical  Bohm 
diffusion  uses  a  =  1 .  HPHALL  was  originally  created  with  a  single  value  for  a  and  then  moved 
to  a  two-region  model,  with  separate  values  of  ac  and  ap  for  the  channel  and  the  plume, 
respectively,  after  it  was  found  that  a  two-region  model  shows  an  improvement  in  accuracy  on 
similar  but  different  codes  {Hagelaar  et  ah,  2003,  Koo  et  al„  2006}.  A  three-region  model  was 
later  implemented  by  Hofer  et  al.  {Mofer  et  al.,  2008}  in  order  to  more  closely  align  cross-field 
mobility  to  the  distribution  of  the  Hall  parameter  within  the  channel,  separating  the  parameter 
into  etc,  etc!  and  ap  for  the  channel,  acceleration  region,  and  plume,  respectively.  A  visual 
representation  of  the  three-region  mobility  model  can  be  seen  in  Figure  3.  The  previous,  two- 
region  model  can  also  be  visualized  by  combining  Regions  I  and  El  of  Figure  3.  As  used  by  Hofer 
et  al.  and  in  order  to  maintain  consistency,  in  the  present  study  values  of  ac  =  0.044,  ap  =  1  are 
used  for  the  two-region  model  and  etc  _  0.08,  ae  =  0.016,  ap  =  10  are  used  for  the  three-region 
model. 


Figure  3:  Experimental  data  of  Hall  parameter  versus  axial 
position  along  with  representation  of  mobility  region. 


The  simulated  discharge  current  traces  of  the  three-region  mobility  model  produce  more 
easily  distinguishable  features  than  the  two-region  model  results,  shown  in  Figure  4  in  which  the 
three-  and  two-region  mobility  models  are  compared  side  by  side.  The  three  most  profound 
phenomena  found  in  the  three  mobility  region  current  traces  are  the  breathing  mode  oscillation,  a 
higher-frequency  mode  called  the  ' Three-region  artifact",  and  the  discharge  current  lag  from  the 
beam  to  the  anode.  Common  to  both  mobility  models  are  the  lag  and  breathing  mode 
oscillations.  In  addition,  as  mentioned  earlier,  the  three-region  model  predicts  a  standard 
deviation  in  the  discharge  oscillations  of  about  8-25\%  of  the  mean  discharge  current  while  the 
two-region  model  predicts  a  steady  deviation  of  approximately  8\%,  These  relative  magnitudes 
can  be  seen  in  the  side  by  side  comparison.  The  spectral  analyses  of  two-region  and  three-region 
simulations  can  be  found  in  Figure  5  for  four  different  simulation  time  steps,  showing  significant 
differences  in  the  results  produced  by  the  two  mobility  models*  The  time  steps  used  in  this 
parametric  study  are  At  —  5e-9  s,  2.5e-8  s,  5e-8  s,  and  10  e~8  s. 

A  major  feature  found  in  the  three-region  current  traces  is  the  high-frequency  mode  which 
has  a  timescale  of  about  1 .5e-5  s,  corresponding  to  a  frequency  range  of  about  60-80-kHz  which 
can  be  identified  as  the  sharp  peaks  in  the  spectra  of  Figure  5.  This  mode  is  completely  absent  in 
the  two-region  results.  In  addition,  higher  frequency  harmonics  due  to  the  high-frequency  mode 
are  not  observed  in  the  larger  time  step  of  At  -  IOe-8  s,  instead  only  displaying  a  single  spectral 
spike.  This  indicates  that  the  oscillation  mechanism  occurs  at  a  timescale  somewhere  in  between 
At  -  5e-8  s  and  At  =  1 0e-8  s. 
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Figure  4:  Simulated  and  beam  current  profiles  using  three 
region  (left)  and  two  region  (right)  models. 
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Figure  5:  Spectral  analysis  for  two  and  three  region  models; 
At  =  5e-9  (top  left),  2,5e-8  (top  right),  5e-8  (bottom  left),  and 
10e-10  (bottom  right). 


C.  Effects  of  Numerical  Parameters 


Despite  the  similarity  of  the  power  spectra  trends  between  experiment  and  simulation  for 
the  6~kW  Hall  thruster,  it  is  found  that  certain  numerical  parameters  greatly  affect  the  simulated 
time-varying  solutions  (as  previously  portrayed  by  Figure  5)  while  the  time-averaged  solutions 
remain  relatively  unchanged.  In  addition,  these  varying  results  are  mobility-model  specific.  Part 
of  the  evaluation  process  of  the  code  is  to  study  any  sensitivity  of  the  solutions  to  the  numerical 
parameters.  Simulation  particle  population,  time  step,  and  the  acceleration  region  (region  II) 
length  scale  of  the  mobility  model  are  varied  in  order  to  characterize  the  different  time-varying 
solutions. 

Spectral  analysis  is  performed  on  the  discharge  current  signal  at  varying  particle 
populations  and  the  results  are  shown  in  Figure  6  in  which  the  three-region  mobility  model  is 
used.  The  numbers  of  both  neutral  and  singly-charged  xenon  are  varied  and  show  no  change  or 
shift  in  spectral  components  at  the  Hall  thruster  oscillation  modes  of  interest.  Low  frequency 
results  cannot  be  explained  but  seem  to  be  weakly  coupled  to  particle  population.  However,  it  has 
been  observed  experimentally  and  speculated  that  one  source  of  such  low  frequency  oscillations, 
often  called  the  "spoke"  mode,  is  coupled  to  density  nonuniformities  and  ionization  processes 
{James  et  al.,  1966}.  The  spoke  mode  is  impossible  to  capture  in  H  PH  ALL  because  the 
phenomenon  is  inherently  azimuthal  while  the  simulation  operates  in  a  radial-axial,  axisymmetric 
manner.  At  these  varying  ion  and  neutral  populations,  tune-averaged  results,  as  well  as  two-  and 
three-region  results,  vary  by  a  negligible  amount. 
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Figure  6:  A  comparison  of  power  spectra,  varying  singly 
charged  (left)  and  neutral  (right)  particle  populations  using 
three  region  mobility  model. 

Spectral  analysis  is  also  performed  on  the  discharge  current  signals  obtained  with 
different  lime  steps  for  both  two-region  and  three-region  simulations  and  the  results  are  shown  in 
Figure  7.  Al  these  different  time  steps,  time-averaged  results  show  significant  changes  in 
discharge  current  and  thrust  (see  fable  2),  one  of  the  first  indications  of  numerical  sensitivity  in 
the  three-region  mobility  model.  Table  2  also  displays  the  very  small  difference  in  time-averaged 
results  for  varying  electron  subcycle.  Decreasing  time  step  shows  a  significant  increase  in  the 


magnitude  of  the  discharge  current  oscillations  at  higher  frequencies  when  operating  with  the 
three-region  model*  This  is  not  observed  in  the  two-region  simulations  w  ith  the  exception  of  the 
smallest  time  step.  At  =  5e-9  s,  though  the  high-frequency  mode  is  still  not  present.  Furthermore, 
both  the  two-region  and  three-region  models  show  similar  trends  in  changing  time-averaged  data 
with  decreasing  time  step,  the  trend  in  fable  2  showing  increasing  magnitudes  until  the  smallest 
time  step.  Tw'o  simulated  discharge  current  oscillations  operating  with  the  three-region  mobility 
model  and  with  an  order  of  magnitude  difference  of  simulation  time  step  are  shown  side-by-side 
in  Figure  8,  displaying  the  magnitude  increase  and  shift.  It  is  also  observed  that  the  oscillations 
around  the  mean  discharge  current  are  not  symmetric;  the  positive  oscillations  are  of  higher 
magnitude  than  the  negative  oscillations.  The  smaller  time  step  simulations  result  in  generally 
higher  magnitude  oscillations  reaching  peak-to-peak  values  on  the  order  of  about  50%  the  mean 
discharge  current  value,  as  mentioned  before. 
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Figure  1:  A  comparison  of  power  spectra,  varying  time  step 
using  two  region  (left)  and  three  region  (right)  models* 


Figure  8:  Simulated  discharge  current,  using  time  steps  of  5e-S 
(left)  and  5c-9  (right)  at  optimal  operating  condition. 
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10  x  10-8  sec 

19.12 

0.390 

Table  2:  Simulated  time  averaged  performance  data  for  three  region  model 

In  order  to  verify  that  the  simulation  time  steps  are  appropriate,  electron  and  ion 
frequencies  are  calculated  {Nicholson,  1983}  using  Equations  4  and  5.  These  frequencies 
correspond  to  ion  timescales  of  approximately  le-9  s  to  5e-7  s,  seen  in  Figure  9,  and  electron 
timescales  which  are  approximately  500  times  smaller.  Despite  this  difference,  it  is  important  to 
note  that  at  the  larger  electron  subcycle,  Af*  =  0.01  At,  solutions  still  show  little  to  no  disparity. 

More  than  an  order  of  magnitude  difference  in  time  scale  occurs  between  the  bulk  of  the 
discharge  within  the  channel  and  the  near-plume  region.  These  results  show  that  the  two  smallest 
time  steps,  5e-9  s  and  2.5e-8  s,  are  small  enough  to  capture  macroscopic  plasma  phenomena 
throughout  the  whole  domain.  The  5e-8  sec  time  step  seems  small  enough  yet  very  close  to  the 
timescales  within  the  bulk  discharge  in  the  channel.  The  IOe-8  s  time  step  is  too  large  to  resolve 
the  detailed  properties  of  the  channel  discharge  (the  bluest  contours  in  Figure  9).  This  reinforces 
the  previous  observation  that  some  of  the  high-frequency  modes  present  in  the  spectral  analyses 
of  the  smaller  time  steps  are  missing  in  the  10e-8  s  studies. 
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Figure  9:  Two  dimensional,  time  averaged  contours  of  ion  timescales. 
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Figure  10:  Spectral  analyses,  three  region  mobility  model,  for  half  (top  left), 
double  (top  right),  and  four  times  (bottom)  the  original  size  of  region  II. 


Spectral  analysis  is  also  performed  on  the  discharge  current  signals  obtained  with 
different  lengths  of  region  II  of  the  three-region  simulations  and  the  results  are  shown  in  Figure 
10.  At  these  different  lengths,  time-averaged  results  show  significant  changes  in  discharge 
current  and  thrust  (see  Table  3).  Firstly,  it  is  important  to  note  that  without  proper  tuning  of  the 
mobility  coefficient,  a,  time-averaged  data  from  differing  mobility  regions  will  rarely  match 


well.  The  purpose  of  changing  mobility  region  lengths  is  primarily  to  observe  any  response  of  the 
high-frequency  mode  so  that  further  analysis  can  be  made  of  its  physical  nature. 

Changing  the  acceleration  region  to  half  its  original  size  (top  left  of  Figure  10)  results  in  a 
loss  of  the  high-frequency  mode,  possibly  because  of  the  inability  to  resolve  the  physics 
introduced  in  this  region.  Doubling  the  size  of  region  II  (top  right  of  Figure  10)  shifts  the  high- 
frequency  modes  approximately  15  kHz  higher.  Increasing  the  size  of  region  II  to  four  times  the 
original  length  (bottom  of  Figure  10)  creates  a  sharper  spike  and  resonances  at  the  breathing 
mode  and  a  broader  bump  centered  at  about  100  kHz.  The  responses  of  both  time-averaged  and 
time-varying  data  with  varying  mobility  region  size  shows  that  the  physical  nature  of  the  high- 
frequency  mode  might  very  well  be  associated  with  the  physics  captured  by  sizing  the  mobility 
regions  to  experimentally  observed  Hall  thruster  data,  as  previously  portrayed  by  Figure  3. 


Region  II  Size 

h  (A) 

Fth  (X) 

Original 

19.98 

0.391 

Half 

22.48 

0.404 

Double 

17.04 

0.351 

x4 

15.86 

0.342 

Table  3:  Simulated,  time  averaged  performance  data  for  three 
region  model,  varying  acceleration  zones. 

2.4  Conclusions 

Numerical  modeling  of  a  6k  W  Hall  thruster  has  been  performed  using  an  established 
hybrid  fluid-PlC  simulation  code.  Comparisons  of  numerical  results  to  measurements  of  time- 
averaged  properties  such  as  thrust  and  discharge  current  showed  very  good  agreement.  However, 
the  same  solutions  showed  relatively  poor  agreement  with  measurements  of  discharge  current 
oscillations.  In  particular,  the  simulations  predicted  higher  magnitude  oscillations  than  those 
observed  experimentally. 

Subsequent  numerical  studies  investigated  the  sensitivity  of  the  solutions  in  the  time 
domain  to  the  numerical  parameters  and  electron  mobility  models  employed.  It  was  found  that 
the  results  were  insensitive  to  the  number  of  particles  employed.  Significant  sensitivity  of  the 
solutions  to  the  time  step  employed  was  found,  particularly  for  the  three  region  mobility  model. 
In  addition,  it  was  found  that  the  length  scales  employed  by  the  three-region  mobility  model 
affect  the  way  in  which  the  high-frequency  mode  is  captured. 

Future  work  will  focus  on  the  development  of  identifying  guidelines  that  allow  solutions 
to  be  computed  that  are  independent  of  the  numerical  parameters  employed.  In  addition, 
attempts  will  be  made  to  modify  the  three-region  electron  mobility  model  to  make  its  solutions 
less  sensitive  to  numerical  parameters.  The  methods  used  to  simulate  the  detailed  plasma  physics 
within  the  discharge  channel  must  eventually  be  improved  to  result  in  a  more  robust  response  to 
varying  numerical  parameters  and  physical  models. 
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